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More than 99% of the mass of the visible universe is made up of protons and 
neutrons. Both particles are much heavier than their quark and gluon con- 
stituents, and the Standard Model of particle physics should explain this dif- 
ference. We present a full ab-initio calculation of the masses of protons, neu- 
trons and other light hadrons, using lattice quantum chromodynamics. Pion 
masses down to 190 mega electronvolts are used to extrapolate to the physi- 
cal point with lattice sizes of approximately four times the inverse pion mass. 
Three lattice spacings are used for a continuum extrapolation. Our results 
completely agree with experimental observations and represent a quantitative 
confirmation of this aspect of the Standard Model with fully controlled uncer- 
tainties. 

The Standard Model of particle physics predicts a cosmological, quantum chromodynamics 
(QCD)-related smooth transition between a high-temparature phase dominated by quarks and 
gluons and a low-temperature phase dominated by hadrons. The very large energy densities at 
the high temperatures of the early universe have essentially disappeared through expansion and 
cooling. Nevertheless, a fraction of this energy is carried today by quarks and gluons, which are 
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confined into protons and neutrons. According to the mass-energy equivalence, E = m ■ , we 
experience this energy as mass. Because more than 99% of the mass of ordinary matter comes 
from protons and neutrons, and in turn about 95% of their mass comes from this confined 
energy, it is of fundamental interest to perform a controlled, ab initio calculation based on QCD 
to determined the hadron masses. 

QCD is a generalized version of quantum electrodynamics (QED) which describes the elec- 
tromagnetic interactions. The Euclidean Lagrangian with gauge coupling g and a quark mass of 
m can be written as C=-l/ {2g^)TrF^^F^^ + ^/j['y^{df, + A^,) +m\'ip, where Ff,^=df,A^-d^Af,+ 
[Afj,,A,y]. In electrodynamics, the gauge potential A^ is a real valued field, whereas in QCD it 
is a 3x3 matrix field. Consequently, the commutator in F^i, vanishes in QED, but not in QCD. 
The ip fields also have an additional "color" index in QCD, which runs from 1 to 3. Differ- 
ent "flavors" of quarks are represented by independent fermionic fields, with possibly different 
masses. In the work presented here, a full calculation of the light hadron spectrum in QCD, 
only three input parameters are required: the light and strange quark masses and the coupling g. 

The action S of QCD is defined as the four-volume integral of C Green's functions are 
averages of products of fields over all field configurations, weighted by the Boltzmann factor 
exp(— S*). A remarkable feature of QCD is asymptotic freedom, which means that for high 
energies (that is, for energies at least 10 to 100 times higher than that of a proton at rest) the 
interaction gets weaker and weaker (T'2), enabling perturbative calculations based on a small 
coupling parameter. Much less is known about the other side, where the coupling gets large, 
and the physics describing the interactions becomes nonperturbative. To explore the predictions 
of QCD in this nonperturbative regime, the most systematic approach is to discretize (3) the 
above Lagrangian on a hypercubic space-time lattice with spacing a, to evaluate its Green's 
functions numerically and to extrapolate the resulting observables to the continuum (a — > 0). A 
convenient way to carry out this discretization is to place the fermionic variables on the sites of 
the lattice, whereas the gauge fields are treated as 3 x 3 matrices connecting these sites. In this 
sense, lattice QCD is a classical four-dimensional statistical physics system. 

Calculations have been performed using the quenched approximation, which assumes that 
the fermion determinant (obtained after integrating over the ^ fields) is independent of the gauge 
field. Although this approach omits the most computationally demanding part of a full QCD 
calculation, a thorough determination of the quenched spectrum took almost 20 years. It was 
shown ^ that the quenched theory agreed with the experimental spectrum to approximately 
10% for typical hadron masses and demonstrated that systematic differences were observed 
between quenched and two flavor QCD beyond that level of precision (4 Jj). 

Including the effects of the light sea quarks has dramatically improved the agreement be- 
tween experiment and lattice QCD results. Five years ago, a collaboration of collaborations © 
produced results for many physical quantities that agreed well with experimental results. Thanks 
to continuous progress since then, lattice QCD calculations can now be performed with light 
sea quarks whose masses are very close to their physical values ( 7) (though in quite small vol- 
umes). Other calculations, which include these sea-quark effects in the light hadron spectrum, 
have also appeared in the literature (|llll[70][22][22l[22][24][I21[Ml)- However, all of these studies 
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have neglected one or more of the ingredients required for a full and controlled calculation. The 
five most important of those are, in the order that they will be addressed below: 
I. The inclusion of the up (n), down (d) and strange (s) quarks in the fermion determinant with 
an exact algorithm and with an action whose universality class is QCD. For the light hadron 
spectrum, the effects of the heavier charm, bottom and top quarks are included in the coupling 
constant and light quark masses. 

n. A complete determination of the masses of the light ground-state, flavor nonsinglet mesons 
and octet and decuplet baryons. Three of these are used to fix the masses of the isospin averaged 
light (mud) and strange (rris) quark masses and the overall scale in physical units. 
III. Large volumes to guarantee small finite-size effects and at least one data point at a signif- 
icantly larger volume to confirm the smallness of these effects. In large volumes, finite-size 
corrections to the spectrum are exponentially small (fJTlIT^ . As a conservative rule of thumb 
M^L>4, with the pion mass and L the lattice size, guarantees that finite-volume errors in 
the spectrum are around or below the percent level (29). Resonances require special care. Their 
finite volume behavior is more involved. The literature provides a conceptually satisfactory 
framework for these effects (|i 91 1201) which should be included in the analysis. 
rV. Controlled interpolations and extrapolations of the results to physical rriud and (or even- 
tually directly simulating at these mass values). Although interpolations to physical m^, cor- 
responding to Af/c— 495 MeV, are straightforward, the extrapolations to the physical value of 
rriud, corresponding to Af^r^lBS MeV, are difficult. They need computationally intensive calcu- 
lations with reaching down to 200 MeV or less. 

V. Controlled extrapolations to the continuum limit, requiring that the calculations be performed 
at no less than three values of the lattice spacing, in order to guarantee that the scaling region is 
reached. 

Our analysis includes all five ingredients listed above, thus providing a calculation of the 
light hadron spectrum with fully controlled systematics as follows. 

/. Owing to the key statement from renormalization group theory that higher-dimension, 
local operators in the action are irrelevant in the continuum limit, there is, in principle, an un- 
limited freedom in choosing a lattice action. There is no consensus regarding which action 
would offer the most cost-effective approach to the continuum limit and to physical rriud- We 
use an action that improves both the gauge and fermionic sectors and heavily suppresses non- 
physical, ultraviolet modes (|29| . We perform a series of 2-1-1 flavor calculations: that is, we 
include degenerate u and d sea quarks and an additional s sea quark. We fix to its approxi- 
mate physical value. To interpolate to the physical value, four of our simulations were repeated 
with a slightly different m^. We vary niud in a range that extends down to :^190 MeV. 

//. QCD does not predict hadron masses in physical units: only dimensionless combinations 
(such as mass ratios) can be calculated. To set the overall physical scale, any dimensionful 
observable can be used. However, practical issues influence this choice. First of all, it should be 
a quantity that can be calculated precisely and whose experimental value is well known. Second, 
it should have a weak dependence on rriud so that its chiral behavior does not interfere with that 
of other observables. Because we are considering spectral quantities here, these two conditions 
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should guide our choice of the particle whose mass will set the scale. Furthermore, the particle 
should not decay under the strong interaction. On the one hand, the larger the strange content 
of the particle, the more precise the mass determination and the weaker the dependence on niud- 
These facts support the use of the Q baryon, the particle with the highest strange content. On 
the other hand, the determination of baryon decuplet masses is usually less precise than those of 
the octet. This observation would suggest that the S baryon is appropriate. Because both the 
and S are reasonable choices, we carry out two analyses, one with Mq (fi set) and one with M= 
(S set). We find that for all three gauge couplings, 6/(7^=3.3, 3.57 and 3.7, both quantities give 
consistent results, namely: a~0.125, 0.085 and 0.065 fm, respectively. To fix the bare quark 
masses, we use the mass ratio pairs Mt,/Mq^,Mk/Mq, or Mt,/M'e,Mk/M'e. We determine the 
masses of the baryon octet (A^, S, A, S) and decuplet (A, S*, S*, Vt) and those members of the 
light pseudoscalar (vr, K) and vector meson (p, K*) octets that do not require the calculation of 
disconnected propagators. Typical effective masses are shown in Figure [IJ 

///. Shifts in hadron masses due to the finite size of the lattice are systematic effects. There 
are two different effects and we took both of them into account. The first type of volume depen- 
dence is related to virtual pion exchange between the different copies of our periodic system and 
it decreases exponentially with M^L. Using Mt,L>^A results in masses which coincide, for all 
practical purposes, with the infinite volume results [see results, for example, for pions (127 1) and 
for baryons (|22l[23l) ). Nevertheless, for one of our simulation points we used several volumes 
and determined the volume dependence which was included as a (negligible) correction at all 
points (|29l) . The second type of volume dependence exists only for resonances. The coupling 
between the resonance state and its decay products leads to a non-trivial level structure in finite 
volume. Based on (79ll20l). we calculated the corrections necessary to reconstruct the resonance 
masses from the finite volume ground-state energy and included them in the analysis (29). 

IV. Though important algorithmic developments have taken place recently [for example ([2?1 
\25\i and for our setup ([231)1, simulating directly at physical m^^ in large enough volumes, which 
would be an obvious choice, is still extremely challenging numerically. Thus, the standard 
strategy consists of performing calculations at a number of larger mud and extrapolating the 
results to the physical point. To that end we use chiral perturbation theory and/or a Taylor 
expansion around any of our mass points ([29l) . 

V. Our three-flavor scaling study ([231) showed that hadron masses deviate from their con- 
tinuum values by less than approximately 1% for lattice spacings up to a!^0.125 fm. Because 
the statistical errors of the hadron masses calculated in the present paper are similar in size, we 
do not expect significant scaling violations here. This is confirmed by Figured Nevertheless, 
we quantified and removed possible discretization errors by a combined analysis using results 
obtained at three lattice spacings ([291 ). 

We performed two separate analyses, setting the scale with M= and Mq. The results of these 
two sets are summarized in Table [TJ The E set is shown in Figure |3l With both scale-setting 
procedures we find that the masses agree with the hadron spectrum observed in nature (27). 

Thus, our study strongly suggests that QCD is the theory of the strong interaction, at low 
energies as well, and furthermore that lattice studies have reached the stage where all systematic 
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errors can be fully controlled. This will prove important in the forthcoming era in which lattice 
calculations will play a vital role in unraveling possible new physics from processes which are 
interlaced with QCD effects. 
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Table 1 : Spectrum results in giga electronvolts. The statistical (SEM) and systematic uncertain- 
ties on the last digits are given in the first and second set of parentheses, respectively. Exper- 
imental masses are isospin-averaged (l29l) . For each of the isospin multiplets considered, this 
average is within at most 3.5 MeV of the masses of all of its members. As expected the octet 
masses are more accurate than the decuplet masses, and the larger the strange content the more 
precise is the result. As a consequence the A mass determination is the least precise. 
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Figure 1 : Effective masses aM=log [C{t/a)/C{t/a + l)], where C{t/a) is the correlator at time 
t, for TT, i^', A^, S and f2 at our lightest simulation point with M^^190 MeV (a ^ 0.085 fm with 
physical strage quark mass). For every 10th trajectory, the hadron correlators were computed 
with Gaussian sources and sinks whose radii are approximately 0.32 fm. The data points rep- 
resent mean ± SEM. The horizontal lines indicate the masses ± SEM obtained by performing 
single mass correlated cosh/sinh fits to the individual hadron correlators with a method similar 
to that of S28i). 
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Figure 2: Pion mass dependence of the nucleon (A^) and Q for all three values of the lattice 
spacing. (A): masses normalized by M=, evaluated at the corresponding simulation points. (B): 
masses in physical units. The scale in this case is set by Ms at the physical point. Triangles 
on dotted lines correspond to a^O.125 fm, squares on dashed lines to a^O.085 fm and circles 
on solid lines to a^O.065 fm. The points were obtained by interpolating the lattice results 
to the physical (defined by setting IMj^-M^ to its physical value). The curves are the 
corresponding fits. The crosses are the continuum extrapolated values in the physical pion mass 
limit. The lattice-spacing dependence of the results is barely significant statistically despite 
the factor of 3.7 separating the squares of the largest (a~0.125 fm) and smallest (a~0.065 fm) 
lattice spacings. The x^/degrees of freedom values of the fits in (A) are 9.46/14 (0,) and 7.10/14 
(A^), whereas those of the fits in (B) are 10.6/14 (Q) and 9.33/14 (A^). All data points represent 
mean ± SEM. 
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Figure 3: The light hadron spectrum of QCD. Horizontal lines and bands are the experimental 
values with their decay widths. Our results are shown by solid circles. Vertical error bars 
represent our combined statistical (SEM) and systematic error estimates, tt, K and S have no 
error bars, because they are used to set the light quark mass, the strange quark mass and the 
overall scale, respectively. 
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Details of the simulations 

We use a tree-level, 0(a^)-improved Symanzik gauge action (57) and work with tree-level, 
clover-improved Wilson fermions, coupled to links which have undergone six levels of stout 
link averaging (S2). (The precise form of the action is presented in {S3).) 

Simulation parameters, lattice sizes and trajectory lengths after thermalization are summa- 
rized in Table ISTl Note, that we work on spatial volumes as large as L^~(4fm)^ and temporal 
extents up to T~8 fm. Besides significantly reducing finite- volume corrections, this choice has 
a similar effect on the statistical uncertainties of the results as increasing the number of trajecto- 
ries at fixed volume. For a given pion mass, this increase is proportional to the ratio of volumes. 
Thus, for T oc L, 1,300 trajectories at M^L=4 are approximately equivalent to 4,000 trajec- 
tories at M^L=3. (A factor comes from the summation over the spatial volume required to 
project the hadron correlation functions onto the zero-momentum sector and an additional L 
comes from the fact that more timeslices are available for extracting the corresponding hadron 
mass.) 

The integrated autocorrelation times of the smeared plaquette and that of the number of 
conjugate gradient iteration steps are less than approximately ten trajectories. Thus every tenth 
trajectory is used in the analysis. We calculate the spectrum by using up to eight timeslices 
as sources for the correlation functions. For the precise form of the hadronic operators see 
e.g. (50). We find that Gaussian sources and sinks of radii ~ 0.32 fm are less contaminated by 
excited states than point sources/sinks (see Figure [STI). The integrated autocorrelation times for 
hadron propagators, computed on every tenth trajectory, are compatible with 0.5 and no further 
correlations were found through binning adjacent configurations. In order to exclude possible 
long-range correlations in our simulations, we performed a run with 10,000 and one with 4,500 
trajectories. No long-range correlations were observed. Further, we never encountered algo- 
rithmic instabilities as illustrated by the time history of the fermionic force in Figure [S2l and 
discussed in more detail in {S3). Note that the fermionic force, which is the derivative of the 
fermionic action with respect to the gauge field, is directly related to the locality properties of 
our action (see Figure IS3l). 

Finite volume corrections and resonances 

For fixed bare parameters (gauge coupling, light quark mass and strange quark mass), the ener- 
gies of the different hadronic states depend on the spatial size of the lattice (in a finite volume 
the energy spectrum is discrete and all states are stable). There are two sources of volume de- 
pendence, which we call type I and type II. These were discussed in a series of papers by M. 
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Liischer ^ Both effects were quantified in a self-consistent manner in our analysis, 

using only the results of our calculations (i.e. no numerical inputs from experiments were used). 

Type I effects result from virtual pion exchanges between the different copies of our peri- 
odic system. These effects induce corrections in the spectrum which fall off exponentially with 
Mt,L for large enough volumes (S5). For one set of parameters (M^^320 MeV at a^O.125 fm), 
additional runs have been carried out for several spatial volumes ranging from M^L^3.5 to 
7. The size dependences of the different hadron masses Mx are successfully described by 
Mx{L) = Mx + cx{M^) ■ exp(-M^L)/(M^L)3/2 Figure |M1 shows the volume dependence 
at M7r=320 MeV for the two statistically most significant channels : the pion and nucleon chan- 
nels. The fitted cx coefficients are in good agreement with those suggested by 6 [70l) which 
predicts a behavior of cx(M^) oc M^. Our results for these and other channels confirm the rule 
of thumb: Mt,L>A gives the infinite volume masses within statistical accuracy. Nevertheless, 
we included these finite volume corrections in our analysis. 

The other source of volume dependence (type II) is relevant only to resonant states, in 
regions of parameter space where they would decay in infinite volume (five out of the twelve 
particles of the present work are resonant states). Since in this case the lowest energy state with 
the quantum numbers of the resonance in infinite volume is a two particle scattering state, we 
need to take the effects of scattering states into account in our analysis. For illustration we start 
by considering the hypothetical case where there is no coupling between the resonance (which 
we will refer to as "heavy state" in this paragraph) and the scattering states. In a finite box 
of size L, the spectrum in the center of mass frame consists of two particle states with energy 

^/Mf+^ + ^Jm^ + k2, where k = n27r/L, n G and Mi, M2 are the masses of the lighter 
particles (with corrections of type I discussed in the previous paragraph) and, in addition, of the 
state of the heavy particle Mx (again with type I corrections). As we increase L, the energy of 
of any one of the two particle states decreases and eventually becomes smaller than the energy 
Mx of X. An analogous phenomenon can occur when we fix L but reduce the quark mass (the 
energy of the two light particles changes more than Mx)- In the presence of interactions, this 
level crossing disappears and, due to the mixing of the heavy state and the scattering state, an 
avoided level crossing phenomenon is observed. Such mass shifts due to avoided level crossing 
can distort the chiral extrapolation of hadron masses to the physical pion mass. 

The literature {!^ ^ ^ provides a conceptually satisfactory basis to study resonances 
in lattice QCD: each measured energy corresponds to a momentum, |k|, which is a solution 
of a complicated non-linear equation. Though the necessary formulae can be found in the 
literature (cf. equations (2.7, 2.10-2.13, 3.4, A3) of (SS)), for completeness the main ingredients 
are summarized here. We follow (S8) where the p-resonance was taken as an example and 
it was pointed out that other resonances can be treated in the same way without additional 
difficulties. The p-resonance decays almost exclusively into two pions. The absolute value 
of the pion momentum is denoted hy k = |k|. The total energy of the scattered particles is 
W = 2(M^ + k'^y^'^ in the center of mass frame. The tttc scattering phase 6ii{k) in the isospin 
/ = 1, spin J = 1 channel passes through 7r/2 at the resonance energy, which correspond to 
a pion momentum k equal to kp = — M^)^/^. In the effective range formula (k^/W) ■ 
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cot^ii = a + bk'^, this behavior implies a = —bk^ = Ak^/{MpTp), where Tp is the decay 
width the resonance (which can be parametrized by an effective coupling between the pions 
and the p). The basic result of {S7) is that the finite- volume energy spectrum is still given by 
W = 2(M^ + k'^y^'^ but with k being a solution of a complicated non-linear equation, which 
involves the nn scattering phase 5u{k) in the isospin 1=1, spin J = 1 channel and reads 
nn — 5ii{k) = (f){q). Here k is in the range < A; < ^/SM^^, n is an integer, q = kL/(2n) and 
(j){q) is a known kinematical function which we evaluate numerically for our analysis (0(g) oc q^ 
for small q and (f)(q) ^ vrg^ for g > 0.1 to a good approximation; more details on 0(g) are given 
in Appendix A of (S8)). Solving the above equation leads to energy levels for different volumes 
and pion masses (for plots of these energy levels, see Figure 2 of {S8)). 

Thus, the spectrum is determined by the box length L, the infinite volume masses of the 
resonance Mx and the two decay products Mi and M2 and one parameter, gx, which describes 
the effective coupling of the resonance to the two decay products and is thus directly related to 
the width of the resonance. In the unstable channels our volumes and masses result in resonance 
states Mx which have lower energies than the scattering states (there are two exceptions, see 
later). In these cases Mx can be accurately reconstructed from L, Mi, M2 and gx- However, 
since we do not want to rely on experimental inputs in our calculations of the hadron masses, 
we choose to use, for each resonance, our set of measurements for various L, Mi and M2 to 
determine both Mx and gx- With our choices of quark masses and volumes we find despite 
limited sensitivity to the resonances' widths, that we can accurately determine the resonances' 
masses. Moreover, the finite volume corrections induced by these effects never exceed a few 
percent. In addition, the widths obtained in the analysis are in agreement with the experimental 
values, albeit with large errors. (For a precise determination of the width, which is not our goal 
here, one would preferably need more than one energy level obtained by cross -correlators. Such 
an analysis is beyond the scope of the present paper.) 

Out of the 14-12=168 mass determinations (14 sets of lattice parameters/volumes-see Ta- 
ble [Sl]-and 12 hadrons) there are two cases for which Mx is larger than the energy of the 
lowest scattering state. These exceptions are the p and A for the lightest pion mass point at 
a~0.085 fm. Calculating the energy levels according to (ST^ ^) for these two isolated cases, 
one observes that the energy of the lowest lying state is already dominated by the contribution 
from the neighboring, two particle state. More precisely, this lowest state depends very weakly 
on the resonance mass, which therefore cannot be extracted reliably. In fact, an extraction of 
Mx from the lowest lying state would require precise information on the width of the reso- 
nance. Since one does not want to include the experimental width as an input in an ab initio 
calculation, this point should not be used to determine Mp and Ma- Thus, for, and only for the 
p and A channels, we left out this point from the analysis. 



13 



Approaching the physical mass point and the continuum Hmit 



We consider two different paths, in bare parameter space, to the physical mass point and contin- 
uum limit. These correspond to two different ways of normalizing the hadron masses obtained 
for a fixed set of bare parameters. For both methods we follow two strategies for the extrap- 
olation to the physical mass point and apply three different cuts on the maximum pion mass. 
We also consider two different parameterizations for the continuum extrapolation. All residual 
extrapolation uncertainties are accounted for in the systematic errors. We carry out this analysis 
both for the E and for the Q sets separately. 

We call the two ways of normalizing the hadron masses: 1. "the ratio method", 2. "mass 
independent scale setting". 

1 . The ratio method is motivated by the fact that in QCD one can calculate only dimen- 
sionless combinations of observables, e.g. mass ratios. Furthermore, in such ratios cancella- 
tions of statistitical uncertainties and systematic effects may occur. The method uses the ratios 
rx=Mx/M^ and parametrizes the mass dependence of these ratios in terms of rj,=MJM'^ and 
rK=MK/M^. The continuum extrapolated two-dimensional surface rx=rx(rTT,i^K) is an unam- 
biguous prediction of QCD for a particle of type X (a couple of points of this surface have been 
determined in {S3)). One-dimensional slices (2r|- — was set to 0.27, to its physical value) of 
the two-dimensional surfaces for N and Q are shown on Figure 2 of our paper. (Here we write 
the formulas relevant for H set; analogous expressions hold for the set. The final results are 
also given for the f2 set). 

A linear term in rj^ (or M|-) is sufficient for the small interpolation needed in the strange 
quark mass direction. On the other hand, our data is accurate enough that some curvature with 
respect to (or M^) is visible in some channels. In order to perform an extrapolation to the 
physical pion mass one needs to use an expansion around some pion mass point. This point can 
be r^=0 (Mt,=0), which corresponds to chiral perturbation theory. Alternatively one can use a 
non-singular point which is in a range of (or M^) which includes the physical and simulated 
pion masses. We follow both strategies (we call them "chiral fit" and "Taylor fit", respectively). 

In addition to a linear expression in M^, chiral perturbation theory predicts (Sll ) an 
next-to-leading order behavior for masses other than those of the pseudo-Goldstone bosons. 
This provides our first strategy ("chiral fit"). A generic expansion of the ratio rx around a 
reference point reads: rx = rx{ref) + ax[rl — r^(re/)] + f3x[rK ~ ''"Ki'i^^f)] + hoc, where 
hoc denotes higher order contributions. In our chiral fit, hoc is of the form r^, all coefficients are 
left free and the reference point is taken to be r^(re/)=0 and r|-(re/) is the midpoint between 
our two values of r^, which straddle r\{j)hys). The second strategy is a Taylor expansion in 
and r\- around a reference point which does not correspond to any sort of singularity ("Taylor 
fit"). In this case, r\{ref) is again at the center of our fit range and r1{rcf) is the midpoint 
of region defined by the physical value of the pion mass and the largest simulated pion mass 
considered. This choice guarantees that all our points are well within the radius of convergence 
of the expansion, since the nearest singularities are at = and/or Mk = 0. Higher order 
contributions, hoc, of the form turned out to be sufficient. 
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We extrapolate to the physical pion mass following both strategies (cubic term of the "chiral 
fit" or a quartic contribution of the "Taylor fit"). The variations in our results which follow from 
the use of these different procedures are included in our systematic error analysis. 

The range of applicability of these expansions is not precisely known a priori. In case 
of the two vector mesons the coefficients of the higher order (r^ or r^) contributions were 
consistent with zero even when using our full pion mass range. Nevertheless, they are included 
in the analysis. For the baryons, however, the higher order contributions are significant. The 
difference between the results obtained with the two approaches gives some indication of the 
possible contributions of yet higher order terms not included in our fits. To quantify these 
contributions further, we consider three different ranges of pion mass. In the first one we include 
all 14 simulation points, in the second one we keep points upto = 0.38 (thus dropping two 
pion mass points) and in the third one we apply an even stricter cut at r^, = 0.31 (which 
corresponds to omitting the five heaviest points). The pion masses which correspond to these 
cuts will be given shortly. The differences between results obtained using these three pion mass 
ranges are included in the systematic error analysis. 

To summarize, the "ratio method" uses the input data rx, r^, and to determine rx{ref ), 
ax and Px and, based on them, we obtain rx at the physical point. The determination of this 
value is done with the two fit strategies ("chiral" and "Taylor") for all three pion mass ranges. 

2. The second, more conventional method ("mass independent scale setting") consists of 
first setting the lattice spacing by extrapolating to the physical point, given by the physical 
ratios of Mt^/M^ and Mk/M^. Using the resulting lattice spacings obtained for each bare gauge 
coupling, we then proceed to fit Mx vs. and applying both extrapolation stratagies 
("chiral" and "Taylor") discussed above. We use the same three pion mass ranges as for the 
"ratio method": in the first all simulation points are kept, in the second we cut at Af^=560 MeV 
and the third case this cut was brought down to M^=450 MeV. 

As shown in the 2+1 flavor scaling study of (^, typical hadron masses, obtained in cal- 
culations which are performed with our 0(a) -improved action, deviate from their continuum 
values by less than approximately 1% for lattice spacings up to a ^ 0.125fm. Moreover, (5(21) 
shows that these cutoff effects are linear in as is scaled from a ~ 0.065 fm to a ~ 0.125 fm 
and even above. Thus, we use the results obtained here, for three values of the lattice spacing 
down to a ~ 0.065 fm, to extrapolate away these small cutoff effects, by allowing rx{ref) 
(or Mx{ref)) to acquire a linear dependence in a^. In addition to the extrapolation in a^, we 
perform an extrapolation in a and use the difference as an estimate for possible contributions of 
higher order terms not accounted for in our continuum extrapolation. 

The physical mass and continuum extrapolations are carried out simultaneously in a com- 
bined, correlated analysis. 



15 



Statistical and systematic error analysis 



Systematic uncertainties are accounted for as described above. In addition, to estimate the possi- 
ble contributions of excited states to our extraction of hadron masses from the time-dependence 
of two-point functions, we consider 18 possible time intervals whose initial time varies from 
low values, where excited states may contribute, to higher values, where the quality of fit clearly 
indicate the absence of such contributions. 

Since the light hadron spectrum is known experimentally it is of extreme importance to 
carry out a blind data analysis. One should avoid any arbitrariness related e.g. to the choice 
of some fitting intervals or pre-specified coefficients of the chiral fit. We follow an extended 
frequentist's method (l ST2\i . To this end we combine several possible sets of fitting procedures 
(without imposing any additional information for the fits) and weight them according to their 
fit quality. Thus, we have 2 normalization methods, 2 strategies to extrapolate to the physical 
pion mass, 3 pion mass ranges, 2 different continuum extrapolations and 18 time intervals for 
the fits of two point functions, which result in 2-2-3-2- 18=432 different results for the mass of 
each hadron. 

In lattice QCD calculations, electromagnetic interactions are absent and isospin is an exact 
symmetry. Electromagnetic and isospin breaking effects are small, typically a fraction of 1% in 
the masses of light vector mesons and baryons {S16). Moreover, electromagnetic effects are a 
small fraction of the mass difference between the members of a same isospin multiplet (j[731). 
We account for these effects by isospin averaging the experimental masses to which we compare 
our results. This eliminates the leading isospin breaking term, leaving behind effects which are 
only a small fraction of 1%. For the pion and kaon masses, we use isospin averaging and 
Dashen's theorem {ST7), which determines the leading order electromagnetic contributions to 
these masses. Higher order corrections, which we neglect in our work, are expected to be below 
the 3 per mil level (see e.g. (.5(7^). All of these residual effects are very small, and it is safe to 
neglect them in comparing our results to experiment. 

The central value and systematic error bar for each hadron mass is determined from the dis- 
tribution of the results obtained from our 432 procedures, each weighted by the corresponding 
fit quality. This distribution for the nucleon is shown in Figure [S5l The central value for each 
hadron mass is chosen to be the median of the corresponding distribution. The systematic error 
is obtained from the central 68% confidence interval. To calculate statistical errors, we repeat 
the construction of these distributions for 2000 bootstrap samples. We then build the bootstrap 
distribution of the medians of these 2000 distributions. The statistical error (SEM) on a hadron 
mass is given by the central 68% confidence interval of the corresponding bootstrap distribu- 
tion. These systematic and statistical errors are added in quadrature, yielding our final error 
bars. The individual components of the total systematic error are given in Table [S2l 
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Table SI: Bare lagrangian parameters, lattice sizes and statistics. The table summarizes the 14 
simulation points at three different lattice spacings ordered by the light quark masses. Note that 
due to the additive mass renormalization, the bare mass parameters can be negative. At each 
lattice spacing 4-5 light quark masses are studied. The results of all these simulations are used 
to perform a combined mass and continuum extrapolation to the physical point. In addition, for 
one set of Lagrangian parameters, different volumes were studied and four of our simulations 
at P=3.57 were repeated with different strange quark masses. 
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Figure SI: Effective masses for different source types in the pion (left panel) and nucleon (right 
panel) channels. Point sources have vanishing extents, whereas Gaussian sources, used on 
Coulomb gauge fixed configurations have radii of approximately 0.32 fm. Clearly, the extended 
sources/sinks result in much smaller excited state contamination. 
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Table S2: Error budget given as fractions of the total systematic error. Results represent av- 
erages over the S and Vt sets. The columns correspond to the uncertainties related to the con- 
tinuum extrapolation {0{a) or 0{a'^) behavior), to the extrapolation to the physical pion mass 
(obtained from chiral/Taylor extrapolations for each of three possible pion mass intervals using 
the ratio method or the mass independent scale setting), to possible excited state contamination 
(obtained from different fit ranges in the mass extractions), and to finite volume corrections 
(obtained by including or not including the leading exponential correction). If combined in 
quadrature, the individual fractions do not add up to exactly 1. The small (<20%) differences 
are due to correlations, the non-Gaussian nature of the distributions and the fact that the very 
small finite volume effects are treated like corrections in our analysis, not contributions to the 
systematic error (the effect of yet higher order corrections is completely negligible). The finite 
volume corrections of the decuplet resonances increase with increasing strange content. This is 
only due to the fact that these are fractions of decreasing total systematic errors. The absolute 
finite volume corrections of these resonances are on the same level. 
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Figure S2: Forces in the molecular dynamics time history. We show here this history for a 
typical sample of trajectories after thermalization. Since the algorithm is more stable for large 
pion masses and spatial sizes, we present -as a worst case scenario- the fermionic force for our 
smallest pion mass (Mtt^IOO MeV; M^L!^4). The gauge force is the smoothest curve. Then, 
from bottom to top there are pseudofermion 1, 2, the strange quark and pseudofermion 3 forces, 
in order of decreasing mass. No sign of instability is observed. 
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Figure S3: Locality properties of the Dirac operator used in our simulations. In the literature, 
the term locality is used in two different ways (see e.g. (ST3\ ST4\ .SiTTI)). Our Dirac operator is 
ultralocal in both senses. First of all (type A locality), in the sum J2x,yi'i^)D{x,y)'ijj{y) the 
non-diagonal elements of our D(x, y) are by definition strictly zero for all (x, y) pairs except 
for nearest neighbors. The figure shows the second aspect of locality (type B), i.e., how D{x, y) 
depends on the gauge field at some distance z: \\dD{x,y) / dU^{x + z)\\. In the analyses 
we use the Euclidian metric for \z\. We take the Frobenius norm of the resulting antihermitian 
matrix and sum over spin, color and Lorentz indices. An overall normalization is performed 
to ensure unity at |^|=0. The action is by definition ultralocal, thus \\dD{x,y) / dU^{x + z)\\ 
depends only on gauge field variables residing within a fixed range. Furthermore, within this 
ultralocality range the decay is, in very good approximation, exponential with an effective mass 
of about 2.2a^. This is much larger than any of our masses, even on the coarsest lattices. 
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Figure S4: Volume dependence of the tt (left panel) and (right panel) masses for one of our 
simulation points corresponding to a ^ 0.125fm and ^ 320 MeV. The results of fits to the 
form C1+C2 exp(— M^L) / (M^L)^/^ are shown as the solid curves, with ci = aMx{L = 00) and 
C2 = acx(M^) given in the text (X = tt,N for pion/nucleon). The dashed curves correspond 
to fits with the cs of refs. 
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Figure S5: Distribution used to estimate the central value and systematic error on the nucleon 
mass. The distribution was obtained from 432 different fitting procedures as explained in the 
text. The median is shown by the arrow. The experimental value of the nucleon mass is indicated 
by the vertical line. 
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